%% Reproducing Results, Figures and Table
% Cleaning the workspace
close all
clear
clc
format short

%% Loading the experimental dataset
load DataSet.mat

% Loading the model data for the four groups
for group=1:4
    for individual=1:ngen{group}
       CR{group,individual}=load(sprintf('SAVED_FILES_Generation_%i_Chromosome_%i/CR.dat',group,individual-1));
    end 
end

%% Computing the number of CRs and CR% for both experimental and model datasets
TrialDuration=1.0; % each trial lasts 1 second
NTrial=77;         % 77 trials per simulation

% For model dataset, for each individual for each group create an array of 77 binary values
% 1 if a CR was identified at that trial, 0 otherwise.
% Data will be stored in CR_MOD cell
for group=1:4
    for individual=1:ngen{group}
        for i=1:NTrial
            % check if a CR was identified at the i_th trial
            greater=CR{group,individual}>=i;
            less=CR{group,individual}>=i-TrialDuration;
            if (sum(less)-sum(greater)==1)
                CR_MOD{group,individual}(i)=1;
            else
                CR_MOD{group,individual}(i)=0;
            end
        end
    end
end

clear CR

% For both model and experimental dataset, compute the CR% along the trials with a mobile
% window of 10 trials
for group=1:4
    CR_PERC_MATRIX_MOD{group}=[];
    for individual=1:ngen{group}
        for i=1:NTrial
            lower=max(1,i-9);
            CR_PERC_MOD{group,individual}(i)=sum(CR_MOD{group,individual}(lower:i))*10; %in percentage
        end
        CR_PERC_MATRIX_MOD{group}=[CR_PERC_MATRIX_MOD{group};CR_PERC_MOD{group,individual}];
    end
end
for group=1:4
    CR_PERC_MATRIX_EXP{group}=[];
    for individual=1:ngen_exp{group}
        for i=1:NTrial
            lower=max(1,i-9);
            CR_PERC_EXP{group,individual}(i)=sum(CR_EXP{group,individual}(lower:i))*10; %in percentage
        end
        CR_PERC_MATRIX_EXP{group}=[CR_PERC_MATRIX_EXP{group};CR_PERC_EXP{group,individual}];
    end
end
clear lower
% For both model and experimental dataset, compute the number (as a percentage) of CRs in EA, LA and EX
for group=1:4
    for individual=1:ngen{group}
        % During EA (Early Acquisition)
        CR_EA_MOD{group}(individual)=sum(CR_MOD{group,individual}(1:11))/11*100;
        % During LA (Late Acquisition)
        CR_LA_MOD{group}(individual)=sum(CR_MOD{group,individual}(57:66))/11*100;
        % During EX (Extinction)
        CR_EX_MOD{group}(individual)=sum(CR_MOD{group,individual}(67:77))/11*100;
    end
end
for group=1:4
    for individual=1:ngen_exp{group}
        % During EA (Early Acquisition)
        CR_EA_EXP{group}(individual)=sum(CR_EXP{group,individual}(1:11))/11*100;
        % During LA (Late Acquisition)
        CR_LA_EXP{group}(individual)=sum(CR_EXP{group,individual}(57:66))/11*100;
        % During EX (Extinction)
        CR_EX_EXP{group}(individual)=sum(CR_EXP{group,individual}(67:77))/11*100;
    end
end

clear greater group i individual less upper

%% Generation of Figure 2a
figure(1)
colors_solid={'k','r','g','b'};
colors_dashed={'--k','--r','--g','--b'};
for group=1:4
    subplot(2,2,group)
    stdshade_median(CR_PERC_MATRIX_MOD{group},0.2,colors_solid{group},[1:NTrial]')
    hold on
    stdshade_median(CR_PERC_MATRIX_EXP{group},0.2,colors_dashed{group},[1:NTrial]')
    set(gca,'fontsize',14);
    box off
    axis([1 NTrial 0 100 ])
    xlabel('Trial','fontsize',14)
    ylabel('%CR','fontsize',14)
end


figure(2)
for group=1:4
    subplot(2,2,group)
    boxplot([CR_EA_EXP{group},CR_EA_MOD{group},...
             CR_LA_EXP{group},CR_LA_MOD{group},...
             CR_EX_EXP{group},CR_EX_MOD{group}], ...
            [zeros(size(CR_EA_EXP{group})),ones(size(CR_EA_MOD{group})),...
             2*ones(size(CR_LA_EXP{group})),3*ones(size(CR_LA_MOD{group})),...
             4*ones(size(CR_EX_EXP{group})),5*ones(size(CR_EX_MOD{group}))], ...
             'Labels',{'EA MOD','EA EXP','LA MOD','LA EXP','EX MOD','EX EXP'},...
             'color',colors_solid{group},'Notch','on');
    set(gca,'XTickLabelRotation',90)
end

%% Computation of Correlation coefficients and the relative p-values
for group=1:4
    [RHO(group,1),RHO(group,2)]=corr(median(CR_PERC_MATRIX_MOD{group})',median(CR_PERC_MATRIX_EXP{group})');
end
disp('     R     p-value');
disp(RHO)

%% Generation of Table 1
for group=1:4
    p(group,1) = ranksum(CR_EA_MOD{group},CR_EA_EXP{group},'alpha',0.01);
    p(group,2) = ranksum(CR_LA_MOD{group},CR_LA_EXP{group},'alpha',0.01);
    p(group,3) = ranksum(CR_EX_MOD{group},CR_EX_EXP{group},'alpha',0.01);
end
disp('     Table 1');
disp(p)

%% Generation of Figure 2b
figure(4)
figure(5)
figure(6)
for group=1:4
    CR_delta{group,1}=CR_EA_MOD{group};
    CR_delta{group,2}=CR_LA_MOD{group}-CR_EA_MOD{group};
    CR_delta{group,3}=CR_EX_MOD{group}-CR_LA_MOD{group};
    for i=1:3 % 0->EA, EA->LA, LA->EX     
        figure(3+i)
        hold on
        med=median(CR_delta{group,i});
        prc25=prctile(CR_delta{group,i},25);
        prc75=prctile(CR_delta{group,i},75);
        plot([0 1+0.02*group],[0 median(CR_delta{group,i})],'color', colors_solid{group},'linewidth',2)
        errorbar(1+0.02*group,med,med-prc25,prc75-med,'color',colors_solid{group},'linewidth',2)
        ylabel('%CR change','fontsize',14)
        box off
        set(gca,'fontsize',16)
        set(gca,'xtick',[])
        clear med prc25 perc75
    end
end

%% Generation of Table 2
% 0->EA
Data=[CR_delta{1,1},CR_delta{2,1},CR_delta{3,1},CR_delta{4,1}];
Groups=[zeros(size(CR_delta{1,1})),ones(size(CR_delta{2,1})),...
             2*ones(size(CR_delta{3,1})),3*ones(size(CR_delta{4,1}))];
[p,tbl,stats] = kruskalwallis(Data,Groups,'off');
stats.gnames={'ALL-PRE';'SHAM-POST';'RIGHT-POST';'LEFT-POST'};
[c_EA,m]=multcompare(stats,'alpha',0.01,'ctype','bonferroni');

% EA->LA
Data=[CR_delta{1,2},CR_delta{2,2},CR_delta{3,2},CR_delta{4,2}];
[p,tbl,stats] = kruskalwallis(Data,Groups,'off');
stats.gnames={'ALL-PRE';'SHAM-POST';'RIGHT-POST';'LEFT-POST'};
[c_LA,m]=multcompare(stats,'alpha',0.01,'ctype','bonferroni');

% LA->EX
Data=[CR_delta{1,3},CR_delta{2,3},CR_delta{3,3},CR_delta{4,3}];
[p,tbl,stats] = kruskalwallis(Data,Groups,'off');
stats.gnames={'ALL-PRE';'SHAM-POST';'RIGHT-POST';'LEFT-POST'};
[c_EX,m]=multcompare(stats,'alpha',0.01,'ctype','bonferroni');

% Draw Table 2
format short e
TABLE2_EA=[NaN, c_EA(1,6), c_EA(2,6), c_EA(3,6)
        c_EA(1,6), NaN, c_EA(4,6), c_EA(5,6)
        c_EA(2,6), c_EA(4,6), NaN, c_EA(6,6)
        c_EA(3,6), c_EA(5,6), c_EA(6,6), NaN];
disp('0->EA'), disp(TABLE2_EA)
TABLE2_LA=[NaN, c_LA(1,6), c_LA(2,6), c_LA(3,6)
        c_LA(1,6), NaN, c_LA(4,6), c_LA(5,6)
        c_LA(2,6), c_LA(4,6), NaN, c_LA(6,6)
        c_LA(3,6), c_LA(5,6), c_LA(6,6), NaN];
disp('EA->LA'),disp(TABLE2_LA)

TABLE2_EX=[NaN, c_EX(1,6), c_EX(2,6), c_EX(3,6)
        c_EX(1,6), NaN, c_EX(4,6), c_EX(5,6)
        c_EX(2,6), c_EX(4,6), NaN, c_EX(6,6)
        c_EX(3,6), c_EX(5,6), c_EX(6,6), NaN];
disp('LA->EX'), disp(TABLE2_EX)

%% Generation of Figure 3
format short
figure
lim=[[-50,50];[-20,20];[-120,120];[-250,250];[-120,120];[-120,120]];
for parameter=1:6
    for group=1:4
        Param{group}=Parameters{group}(:,parameter);
    end
    MeanParam{1}(parameter)=mean(Param{1});
     % Computed as percentage of the ALL-PRE group mean
    for group=1:4
        Param_perc{group}(parameter,:)=100.*(Param{group}-MeanParam{1}(parameter))./MeanParam{1}(parameter);
    end
    subplot(3,2,parameter)
    hold on
    plot([0 4],[mean(Param_perc{1}(parameter,:)) mean(Param_perc{1}(parameter,:))],'k','linewidth',2);
    plot([0 4],[mean(Param_perc{1}(parameter,:))+std(Param_perc{1}(parameter,:)) mean(Param_perc{1}(parameter,:))+std(Param_perc{1}(parameter,:))],'k','linewidth',1);
    plot([0 4],[mean(Param_perc{1}(parameter,:))-std(Param_perc{1}(parameter,:)) mean(Param_perc{1}(parameter,:))-std(Param_perc{1}(parameter,:))],'k','linewidth',1);
    errorbar(1,mean(Param_perc{2}(parameter,:)),std(Param_perc{2}(parameter,:)),'r.');
    errorbar(2,mean(Param_perc{3}(parameter,:)),std(Param_perc{3}(parameter,:)),'g.');
    errorbar(3,mean(Param_perc{4}(parameter,:)),std(Param_perc{4}(parameter,:)),'b.');
    ylim([lim(parameter,1),lim(parameter,2)]);
    set(gca,'fontsize',16)
end
clear lim

%% Statistical differences between parameters of ALL-PRE, SHAM-POST, RIGHT-POST and LEFT-POST
% (The asterisks of Figure 3(a))
for parameter=1:6
    figure
    Data=[Param_perc{1}(parameter,:),Param_perc{2}(parameter,:),...
          Param_perc{3}(parameter,:),Param_perc{4}(parameter,:)];
    Groups=[zeros(size(Param_perc{1}(parameter,:))),ones(size(Param_perc{2}(parameter,:))),...
             2*ones(size(Param_perc{3}(parameter,:))),3*ones(size(Param_perc{4}(parameter,:)))];
    [p,tbl,stats] = kruskalwallis(Data,Groups,'off');
    stats.gnames={'ALL-PRE';'SHAM-POST';'RIGHT-POST';'LEFT-POST'};
    [c{parameter},m]=multcompare(stats,'alpha',0.01,'ctype','bonferroni');
end
clear c m p tbl stats

%% Statistical differences between parameters of SHAM-POST, RIGHT-POST and LEFT-POST
% (The asterisks of Figure 3(b))
for parameter=1:6
    figure
    Data=[Param_perc{2}(parameter,:),Param_perc{3}(parameter,:),Param_perc{4}(parameter,:)];
    Groups=[zeros(size(Param_perc{2}(parameter,:))),ones(size(Param_perc{3}(parameter,:))),...
             2*ones(size(Param_perc{4}(parameter,:)))];
    [p,tbl,stats] = kruskalwallis(Data,Groups,'off');
    stats.gnames={'SHAM-POST';'RIGHT-POST';'LEFT-POST'};
    [c{parameter},m]=multcompare(stats,'alpha',0.01,'ctype','bonferroni');
end
clear c m p tbl stats

%% Generation of Figure 4
% Histogram of changes of each synaptic weight, between all consecutive
% trials during session-2 in SHAM-POST, for each plasticity site

% Loading Synaptic Weights of all groups
% The loading may take some minutes
InitialWeights{1}=load('WeightsEBCC_0_ALLPRE.dat');
InitialWeights{2}=load('WeightsEBCC_0.dat');
InitialWeights{3}=load('WeightsEBCC_0.dat');
InitialWeights{4}=load('WeightsEBCC_0.dat');

for group=1:4
    for t=1:NTrial
        FinalWeights{group,t}=load(sprintf('SAVED_FILES_Generation_%i_Chromosome_0/FinalWeights%i.dat',900+group,t));
    end
end

% Compute the total number of synapses
total = 0;
for w=1:size(InitialWeights{1},1)
    total=total+InitialWeights{1}(w,1);
end

for group=1:4   
    weights=zeros(total,NTrial+1);
    count = 1;
    
    M = InitialWeights{group};
    for w=1:size(M,1)
        if M(w,1)==1
               weights(count,1)=M(w,2);
               count = count+1;
        else
            while(M(w,1)>=1)
                    weights(count,1)=M(w,2);
                    M(w,1) = M(w,1)-1;
                    count = count+1;
            end
        end
    end
    
    k=2;
    for t=1:NTrial
        count = 1;    
        M = FinalWeights{group,t};
        for w=1:size(M,1)
            if M(w,1)==1
               weights(count,k)=M(w,2);
               count = count+1;
            else
                while(M(w,1)>=1)
                    weights(count,k)=M(w,2);
                    M(w,1) = M(w,1)-1;
                    count = count+1;
                end
            end
        end
        k=k+1;
    end
    
    % Selecting the weights of the three plasticities
    First{group}=weights(24001:196805,:);
    Second{group}=weights(369517:374917,:);
    Third{group}=weights(380317:380352,:);
    clear t weights count w k M
    disp(['Loading completed for group ' num2str(group)])
end

% Compute the changes between two consecutive trials
Delta_First = First{2}(:,2:78)-First{2}(:,1:77);
Delta_First = Delta_First(:);

Delta_Second = Second{2}(:,2:78)-Second{2}(:,1:77);
Delta_Second = Delta_Second(:);

Delta_Third = Third{2}(:,2:78)-Third{2}(:,1:77);
Delta_Third = Delta_Third(:);


% Generate the histograms
figure
subplot(1,3,1)
histogram(Delta_First,10,'facecolor','k')
set(gca,'YScale','log')
xlim([-2,2])
box off
ylabel('# Synaptic changes','fontsize',14)
xlabel('Weight Change [ns]','fontsize',14)
set(gca,'fontsize',14)

subplot(1,3,2)
histogram(Delta_Second,10,'facecolor','k')
set(gca,'YScale','log')
xlim([-12e-6 12e-6])
box off
xlabel('Weight Change [ns]','fontsize',14)
set(gca,'fontsize',14)

subplot(1,3,3)
histogram(Delta_Third,10,'facecolor','k')
set(gca,'YScale','log')
xlim([-12e-6 12e-6])
box off
xlabel('Weight Change [ns]','fontsize',14)
set(gca,'fontsize',14)
%% Generation of Figure 5

% Compute the change of synaptic weights from the i_th trial with respect to the initialization (i=1)
% For PF-PC
figure
hold on
for group = 2:4
    for trial=1:77
        Change_weight_First{group}(:,trial)=sum(First{group}(:,trial+1)-First{group}(:,1));
    end
    BlockChange_First{group}(1)=0.0;
    conta=2;
    for i=1:11:77
        A=Change_weight_First{group}(:,i:i+10);
        A=A(:);
        BlockChange_First{group}(conta)=sum(A(:));
        conta=conta+1;
    end
    plot(BlockChange_First{group},colors_solid{group},'marker','x','linewidth',2)
    box off
    set(gca,'fontsize',20)
    xlim([1 8])
    m_First(group-1,:)=diff(BlockChange_First{group});
end

% For MF-DCN
figure
hold on
for group = 2:4
    for trial=1:77
        Change_weight_Second{group}(:,trial)=sum(Second{group}(:,trial+1)-Second{group}(:,1));
    end
    BlockChange_Second{group}(1)=0.0;
    conta=2;
    for i=1:11:77
        A=Change_weight_Second{group}(:,i:i+10);
        A=A(:);
        BlockChange_Second{group}(conta)=sum(A(:));
        conta=conta+1;
    end
    plot(BlockChange_Second{group},colors_solid{group},'marker','x','linewidth',2)
    box off
    set(gca,'fontsize',20)
    xlim([1 8])
    m_Second(group-1,:)=diff(BlockChange_Second{group});
end

% For PF-PC
figure
hold on
for group = 2:4
    for trial=1:77
        Change_weight_Third{group}(:,trial)=sum(Third{group}(:,trial+1)-Third{group}(:,1));
    end
    BlockChange_Third{group}(1)=0.0;
    conta=2;
    for i=1:11:77
        A=Change_weight_Third{group}(:,i:i+10);
        A=A(:);
        BlockChange_Third{group}(conta)=sum(A(:));
        conta=conta+1;
    end
    plot(BlockChange_Third{group},colors_solid{group},'marker','x','linewidth',2)
    box off
    set(gca,'fontsize',20)
    xlim([1 8])
    m_Third(group-1,:)=diff(BlockChange_Third{group});
end

%% Generating Table 3
disp('    PF-PC')
disp('    m1    m2    m3    m4    m5    m6    m7')
disp(m_First)

disp('    MF-DCN')
disp('    m1    m2    m3    m4    m5    m6    m7')
disp(m_Second)

disp('    PC-DCN')
disp('    m1    m2    m3    m4    m5    m6    m7')
disp(m_Third)

